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ABSTRACT 

The Andromeda galaxy (M31) shows many tidal features in its halo, including the Giant 
Southern Stream (GSS) and a sharp ledge in surface density on its western side (the W Shelf). 
Using DEIMOS on the Keck telescope, we obtain radial velocities of M3 1 's giant stars along 
its NW minor axis, in a radial range covering the W Shelf and extending beyond its edge. 
In the space of velocity versus radius, the sample shows the wedge pattern expected from a 
radial shell, which is detected clearly here for the first time. This confirms predictions from 
an earlier model of formation of the GSS, which proposed that the W Shelf is a shell from 
the third orbital wrap of the same tidal debris stream that produces the GSS, with the main 
body of the progenitor lying in the second wrap. We calculate the distortions in the shelf 
wedge pattern expected from its outward expansion and angular momentum, and show that 
these effects are echoed in the data. In addition, a hot, relatively smooth spheroid population 
is clearly present. We construct a bulge-disk-halo A^-body model that agrees with surface 
brightness and kinematic constraints, and combine it with a simulation of the GSS. From the 
contrasting kinematic signatures of the hot spheroid and shelf components, we decompose the 
observed stellar metallicity distribution into contributions from each component using a non- 
parametric mixture model. The shelf component's metallicity distribution matches previous 
observations of the GSS superbly, further strengthening the evidence they are connected and 
bolstering the case for a massive progenitor of this stream. 

Key words: galaxies: M31 - galaxies: interactions - galaxies: kinematics and dynamics 



1 INTRODUCTION 

The Andromeda galaxy (M31) lies close enough for us to resolve 
its individual stars, at a position relatively unobscured by the Milky 
Way (MW) disk. Thus it represents perhaps our best opportu- 
nity to examine the complete structure of a large spiral galaxy in 
great detail. The region that we might term M31's in ner halo or 
outer disk appears quite irregular in resolved-sta r maps (llbata et al.l 
20011: iFerguson et all |2002|- llrwin et all 120051 ; llbata et al.l l2007t 
McConnachie et alj|2009l ; iTanaka et alj|20ld) . These maps show 
numerous substructures, some of which have been catalogued 
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and named, although stars are distributed continuously throughout 
M31's halo out to at least 150 kpc. 

Pencil-beam spectrosc opy and ph otometry brings some order 
to this complicated scene, llbata et al.l d2005f) argues much of the 
structure is an outer extension of M31's disk, with a clumpy struc- 
ture and hot or thick disk kinematics. It also seems much of the 
struct u re is due to M3 1 's gian t southern stream or GSS (llbata et al.l 
l200lh . IFerguson et al.l d2002h first suggested a connection between 
the GSS and the "NE Shelf" feature extending from the NE side 
of the M31 disk (see Fig.fj}, which appears to have a red popu- 
lation of re d giant bran ch (RGB) stars and thus is inferred to be 
metal-rich. Far dal et ail d2007l hereafter F07) pointed out a sim- 
ilar shelf on the western side of M31 (W Shelf), which they in- 
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f erred to be a third wrap of the stream us i ng an jV-body simulation 
of the stream's formation. Gilb ert et all d2007l) found a structure 
through RGB star kinematics on the southeast minor axis (the "SE 
Shelf") and ident i fied it as the fourth wrap of this stream. Indeed, 
Richardson et al. I d2008t> found most of the dense structures in the 
inner halo examined by their deep Hubble Space Telescope (HST) 
observations could be classified as either "disk-like" or "stream- 
like", with the positions of the "stream-like" fields matching well 
to the F07 model. 

These two components cannot be the whole story, how- 
ever. Spectroscopy also offers evidence of a smoother hot halo 



component dGuhathakurta et al.ll2005ll2006l : lchapman et ai1l2006l ; 
iKalirai et alj 2006al ; iGilbert et al. 2007 ). There appears to be an 



upward break in t he surface brightness profile around 20- 30 kpc 
jlrwin et al.ll2005l ; IGuhathakurta et al.ll2005l : llbata et alj|2007l) and 
a decrease in the metallic ity in about the same place dKalirai et al.1 
l2006al : iKoch et alj 120081) . It is much disputed whether the "in- 
ner hal o" region from 10-20kpc is dominated by a halo c om- 
ponent jPritchet & van den Berghl Il994 iDurrell et al 1200 ll) the 
bulge (IGuhathakurta et al.1 12005 1 : IKalirai et alj l2006aT) . the disk 
dlbata et alJl2007l ; ICourteau et al.1 201 lh . or by a mix of these plus 
cold substructures ( Gil bert et al .120071) . Other kinematic anomalies 
such as the "second stream" detected by IKalirai et alj d2006bl) and 
IGilbert et al.1 d2009h currently lack an explanation. One might also 
expect tidal features from M32 or NGC 205 (M 1 1 0), due to their lo - 
cation near M3 1 's center and distorted isophotes dChoi et alj2002h . 

To better understand this complex structure, it is important to 
deduce how many components there are and determine for each 
its distribution in space, velocity, stellar age, and metallicity. A 
key ingredient in this will be a firm determination of the orbit of 
the GSS and the spatial coverage of its debris. The F07 model 
matches many features of the morphology and kinem atics, but it 



is not the only suggestion for the orbit of the stream ( Ibata et all 
l200ll ; iMerrett et al.fl2003l: iFardal et al.llioOol; lHammer et alj|2010l) . 
A crucial test of the model comes in the W shelf region. The shelf 
stars are expected to exhibit a distinctive kinematic pattern in this 
model. F07 presented some kinematic evidence for this model from 
planetary nebulae and RGB stars, but this was more suggestive than 
conclusive. 

In this paper we present new spectroscopy of the resolved stel- 
lar population in the W Shelf region that shows clear evidence for 
distinct shelf and hot spheroid components. This kinematic data 
agrees very well with the prediction from an updated version of the 
F07 model, which was run before the observations were reduced. 
§2 provides details on the observations and our N-body model. §3 
analyzes the observations, first showing (§3.1) the distributions of 
the stars in position, velocity, and metallicity space. We then use 
the iV-body model (§3.2) to separate the shelf and hot halo compo- 
nents in the observations, and thereby infer their distinct metallicity 
distributions. §4 concludes the paper by connecting our results to 
previous results on M31's GSS and other inner halo populations, 
and placing our results in the broader context of halo structure in 
other galaxies. An Appendix provides a discussion of shell kine- 
matics, including a formula for the velocity envelope in the generi- 
cally correct case of an expanding shell. 




Figure 1. Top: star-count map from the INT survey of M31 dlrwin et all 
image created by M. Irwin), with our Keck/DEIMOS spectroscopic 
multislit masks indicated in red. Blue, green, and red channels indicate the 
density of stars in separate color ranges on the RGB. The red channel con- 
tains the reddest, most metal-rich RGB stars. M3 1 's disk is apparent in red, 
and the main tidal features discussed in the text are indicated by labels. The 
W Shelf in particular is just visible as a faint but sharp edge to the NW of 
M3 1 . The box shows the area covered in the bottom panel. The looping ar- 
row indicates the progenitor's path in our model. Bottom: enlargement of 
the region around our spectroscopic masks. Here a matched-filt er map from 
the Subaru/SuprimeCam imaging survey dTanaka et al.[|2010T) is overlaid, 
forming a diagonal stripe from lower left to upper right. The radial zones 
1-5 used later are indicated. 



2 METHODS 

2.1 Observational data 

The data in this paper was obtained as part of the Spectroscopic and 
Photometric Landscape of Andromeda's Stellar Halo (SPLASH) 



survey. Our spect roscopic survey are a is embedded within the 
imaging survey of iTanaka et al.l d2010h . conducted with Suprime- 
Cam on the Subaru Telescope. That survey had the goal of quanti- 
fying the surface brightness profile and metallicity of the M31 halo 
and detecting streams and other substructure, by means of deep 
color-magnitude diagrams (CMDs) obtained mostly along the SE 
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Table 1. Spectroscopic observations 





Radial 


Radius 


Pointing 


center: 


PA 


Observation 


Exposure 


Seeing 


Science 


Confirmed 


M31 


Mask 


Zone 


(kpc) 


a J2000 


5.12000 


o 


Date (UT) 


Time (s) 


(") 


Targets" 


Stars* 


Stars' 1 


NW1V.1 


1 


13.1 


0:38:36.33 


41:50:11.8 


127 


2009 Aug 22 


3600 


0.5-0.8 


279 


56 


46 


NWldV 


2 


18.1 


0:37:07.94 


42:05:35.8 


-143 


2009 Aug 21 


3600 


0.6-0.7 


213 


102 


78 


NW2V.1 


3 


20.6 


0:36:20.72 


42:12:08.5 


-53 


2009 Aug 21 


3000 


0.5-0.6 


213 


89 


62 


NW2V_2 


3 


21.6 


0:36:07.90 


42:15:48.4 


127 


2009 Aug 21 


3000 


0.5-0.7 


208 


100 


62 


NW2dV 


4 


24.2 


0:35:18.87 


42:23:22.8 


-143 


2009 Aug 22 


3300 


0.5-0.6 


201 


124 


90 


NW3V.1 


5 


26.8 


0:34:22.67 


42:28:13.8 


-53 


2009 Aug 21 


3000 


0.5 


193 


76 


25 



"Total number of slits on the mask. 

^Objects with secure radial velocities identified as either MW dwarfs or M31 giants. 
' Objects with secure radial velocities identified as M31 giants. 



and NW minor axes of M31. The SE minor axis had earlier been 
studied intensively by means of wide-field photometry of individ- 
ual stars, extremely deep CMDs obtained with HST, and resolved- 
star spectroscopy jP ritchet & van den Berghl Il994t I Durrell et al.l 
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2007i 120081 : llbata et al.ll2007l : iGilbert et alj|2007l) . In part because 
of increasing MW foreground to the north, the NW minor axis had 
received much less attention. Along the NW minor axis, the Subaru 
survey consists of a series of partially overlapping images that form 
a continuous stripe in the range llkpc <A< 104 kpc. The survey 
used V and / filters on the Johnson-Cousins Vega-based magnitude 
system and probed down to 50% completeness limits of V = 25.7, 
/ = 24.9. 

Our spectroscopic survey leverages the earlier imaging survey 
to obtain kinematics of M31's inner spheroid region, in particular 
covering the W Shelf feature that is prominent in Fig. Q] The blue, 
green, and red channels of the Subaru map are constructed from 
explicit cuts in metallicity with boundaries in [Fe/H] of —2.31, 
— 1.71, —0.71, and 0.0, and due to their different input data and 
construction do not match the larger-scale Isaac Newton Telescope 
(INT) map exactly. Nevertheless, both maps show the W Shelf 
clearly as a red, metal-rich feature with a sharp edge, which is more 
easily seen in the narrow Subaru portion in the lower panel. The ob- 
servation pattern of the Subaru survey involved small regions where 
adjacent Suprime-Cam pointings fields overlapped, and the result- 
ing source catalog is broken up into regions with either single or 
double coverage. We orient the DEIMOS masks parallel to M31's 
minor axis, except within the overlap regions, where we orient the 
masks perpendicular to the minor axis to fit inside the double cover- 
age area. The six spectroscopic masks are described in Table Q]and 
shown in Fig.Q] Because two of the masks nearly overlap, these 
six masks result in five separated radial "zones", labeled 1-5 from 
innermost to outermost in the discussion below. Comparing to the 
maps from the INT and Subaru imaging surveys shown in Fig. [T] 
we can see that zones 1-4 lie within the boundary of the W Shelf 
feature and zone 5 just outside of it. 

The spectroscopic targets were chosen from th e catalog of 
sourc es obtained from the Subaru survey of M31 (ITanaka et al.l 
2010). This catalog was produced using PSF-fitting photometry us- 
ing DAOPHOT-II. Extended objects (mostly galaxies, as well as 
some blended stars) were removed from the target sample based on 
a set of DAOP HOT-II morphologi cal parameters. Details of this 
procedure are in lTanaka et al.l J2010h - 

At these relatively high surface brightness levels in An- 
dromeda's halo, we expect from previous experience that most of 
the stellar targets are M3 1 stars, at least below the tip of the RGB 
at M31's distance at / = 20.5. For simplicity, then, our target se- 



lection uses only the / magnitude, which roughly measures the flux 
in the wavelength range of our spectra. We choose targets in the 
range 20.5 < / < 23.0. We assign a priority to each target which de- 
creases linearly with magnitude by a factor of two over this range. 
Our mask design software then chooses a slit pattern with the aim 
of maximizing the sum of priorities over the mask using a mini- 
mum slit length criterion. This produces a bias toward selection of 
the brighter stars, though stars are selected in significant numbers 
within the entire magnitude range specified. 

Observations were obtained with DEIMOS at the Keck II tele- 
scope in the fall of 2009 in goo d seeing conditions . The observa- 
tional setup was similar to that in lGilbert et al.ld2007t) . We observed 
3 sub-exposures for each of the six masks, with a total observing 
time of 45 to 60 minutes. Details of the observations are listed in 
Table [T] We used the DEIMOS spectrograph with a 1200mm -1 
grating. The central wavelength was 7800 A, with a range of about 
6450-9150 A (depending somewhat on the slit location on the 
mask). The dispersion was 0.33 A pixel -1 , and the spectral res- 
olution about 1.3 A FWHM. 

The spectra were reduced and analyzed using a modified ver- 
sion of the SPEC2D and SPEC ID software developed by the 
DEEP2 team at the University of California, BerkelejQ. These 
routines perform standard spectral reduction steps, including flat- 
fielding, night-sky emission line removal, and extraction of one- 
dimensional spectra from the two-dimensional spectra. Reduced 
one-dimensional spectra were cross-correlated with a library of 
template stellar spectra to determine the radial velocity of the ob- 
ject. Each spectrum was visually inspected and assigned a quality 
code based on the number and quality of absoiption lines. Spectra 
with at least two spectral features (even if one of them is marginal) 
were considered to have secure velocity measurements. The posi- 
tion of the atmospheric A-band in the observed spectrum was used 
to correct t he observed velocity for imperfect cent ering of the star 
in the slit dSimon & Geha1l2007l : ISohn et alj|2007l) . A heliocentric 
correction was applied to the measured velocities based on the sky 
position of each mask and the date and time of the observation. 
More d etails of these procedures a re p rovided in earlier pap ers, in- 
cluding iGuhathakurta etail ||2006) and -Gilbert et al.l j2007l) . 

On several of the masks, a significant fraction of the slits failed 
to yield acceptable spectra. Most of these failures were localized 
by position, so that broad swaths of the mask area are missing any 
velocity measurements. These failures were caused by temporary 
problems with the CCD readout electronics and possible problems 
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with the calibration data. We are confident that the rate of such fail- 
ures has no intrinsic dependence on velocity, and thus introduces 
no bias into our kinematic analysis below. 

The median velocity error for the stars presented in this paper 
is 3kms _1 , and the maximum is 13kms _1 . Errors are estimated 
from the cross-correlation profile, including a modest scaling factor 
that is empirically calibrated by previous duplicate measurements 
of individual stars on overlapping DEIMOS masks. We ignore these 
small velocity errors in the analysis below, as the intrinsic widths 
of the velocity distributions we are interested in are much larger. 

For our metallicity estimates we use photometric rather than 
spectroscopic methods, since for our combination of deep photom- 
etry and moderate S/N spectroscopy we judge the former to be 
more reliable. We estimate metalliciti es from the V and / magni- 
tudes using the procedure descr ibed in Gil bert et al.l (120071) (rather 
than that in lTanakaetal1l2010h . for consistency with the dwarf- 
giant separation algorithm. The [Fe/H] values are based on a com- 
parison of the star's position within the CMD to a finely spaced grid 
of theoretical 12 Gyr, [ct/Fe] = stellar isochrones ( IKalirai et al.l 
l2006ak IVandenBerg et all 120061) adjusted to the distance of M31. 
Metallicities of the MW stars as assigned by this algorithm are 
therefore not meaningful. 

We classify the st ars in the sample using the likelihood method 
of lGilbert et al. (2006). We use the following criteria: 

- The measured equivalent width of the Na I doublet at 8190 A 
combined with the V — I color of the star. 

- The strength of the Ca II triplet absorption lines, as measured 
by comparison of the star's photometric (CMD-based) vs. spectro- 
scopic (Ca triplet-based) metallicity estimates. 

- The position of the star in an (V,/) CMD with respect to theo- 
retical RGB isochrones. 

- The radial velocity of the star. 

Gilbert et al. showed that this combination of diagnostics can 
yield a rea sonably reliable se paration of MW from M3 1 stars. The 
method of Gilbe rt et al.l d2006h can also incorporate additional cri- 
teria such as photometry in the surface-gravity-sensitive DD051 
band, but this was not available for our sample. We flag stars with 
a M3 1 membership probability of greater than 50% and CMD lo- 
cations reasonably close to the RGB as M31 stars (classes 1-3 in 
the terminology of lGilbert et al 1 l2006l) . The membership probabil- 
ity is quite bimodal, and over the sample of identified M31 stars, 
the mean membership probability is 94%. A few misclassifications 
are inevitable, particularly in the higher-velocity part of the sample 
(v ~ — lOO kms -1 ) where the M 31 and MW velocity distributions 
overlap (see Gilbert et al. 2007 for a full discussion). Our final sam- 
ple contains 547 stars, of which 363 are flagged as M3 1 giants and 
184 as MW dwarfs. 

2.2 Simulation 

The A'-body simulation we use in this paper is a slightly updated 
version of the simulation in F07. From our ongoing work with au- 
tomated fitting of the GSS (Fardal et al., in preparation), we se- 
lected a set of parameters that offered good agreement with a set 
of observational constraints . The orbital initi al conditions in the 
sky coordinate system (see Far dal et al. 1 120061) are x = — 24.0kpc, 
y = 0.09 kpc, z = 44. 1 kpc, v x ■ = 29. 1 kms" 1 , v y = 8.7km s~ 1 , v z = 
— 77.6kms~', and the model was run for 0.972Gyr to reach its 
present point. The input satellite is a Plummer model with mass 
M s = 3.09 x 10 9 Mq and scale length r s = 1.16kpc, created with 
the ZENO library. We take the mass to be purely stellar here. We 



reason that much of the dark halo mass will have been stripped 
off in earlier pericent ric passages, as is ty pically found in cos- 
mological simulatio ns dLibeskind et al.ll201 lh and analytic models 
(Watso neTal.l2012h of satellite accretion. The initial mass indicates 
a relat ively massive galax y rather than a dark-matter-dominated 
dwarf. Mori & Rich (2008) have considered a GSS progenitor with 
a massive halo, finding that a total satellite mass of <5 x 10 9 M Q 
is required to avoid perturbing M31's disk too much. (The effects 
of dark matter in the progenitor will be considered in more depth 
in Fardal et al., in preparation, where we find that the moderate 
amounts of dark matter allowed by various constraints do not radi- 
cally change the nature of the simulations.) 

The progenitor initially follows the trajectory shown in Fig.Q] 
(which assumes a test-particle orbit), and though its tidal debris is 
significantly extended it generally follows a similar trajectory. This 
makes the GSS the first wrap of the orbit since the initial pericenter, 
the NE shelf the second, and the W Shelf the third. The W Shelf 
grows outwards over the last 300-400 Myr from M3 1 's center to its 
current size. For convenience we refer throughout to all the parti- 
cles from this simulation as "GSS" material, except when explicitly 
distinguishing different radial wraps. 

We base our fixed gravitational potential model on 
lGeehanetal] {2006). We use a Hernquist bulge with r\, = 0.61 kpc 
and Mi, = 3.39 x 10 10 Mq, a Miyamoto-Nagai disk (replacing 
our earlier exponential disk for computational speed) with a d = 
5.94kpc, b d = 0.4kpc, and M d = 8.07 x 10 10 M n kp c~ 2 , and a 
Navarro-Frenk- White (NFW) halo dNavarro et"ai1ll996l) with r h = 
24.2kpc, p h = 8 c p c = 5.85 x 10 6 M Q kpc" 3 , and M(r < r h ) = 
2.02 x 10 Mq. Our model is qualitatively quite similar to that 
in F07, but better fits the velocity in the GSS itself and the radii 
of the NE and W shelves. The run was conducted with the code 
PKDGRAV, in advance of the observational reductions^ 

We supplement this GSS simulation with an A'-body realiza- 
tion of M31 itself. We draw parameters for the bulge from the 
model just listed. However, for better agreement with the outer sur- 
face brightness profile of M31, we truncate the bulge density start- 
ing at 11 kpc using the default ZENO truncation scheme, which 
assumes continuous values and slopes at the truncation point and 
thereafter adopts a form p(r) oc r~ 2 exp(— r/r t ). For the disk we 
use a true exponential disk as in the iGeehan et al. I J2006h model, 
with M c i — 7.34 x 10 10 MQkpc~ 2 and r d = 5.4kpc, which is sim- 
ilar to the Miyamoto disk used in the simulation in terms of its 
dyna mical effects. We add a stellar halo component not present 
in the Gee han et al.l |2006) model, postulating that it follows the 
dark halo component and thus has r;, = 24.2kpc. The mass is 
normalized to match surface brightness measurements and kine- 
matic samples as described below. We set the characteristic density 
p = 8 c p c = 1.30 x lO 4 M kpc" 3 ,orM(r<r/,) =4.50x 1O 8 M . 
As with the input satellite, ZENO generates the spheroid compo- 
nent particle velocities from the full distribution function generated 
by inversion of an Abel integral, to ensure equilibrium with the po- 
tential of the baryonic components and dark halo. Although we do 
not evolve the A'-body snapshot, from previous test simulations we 
know it is very close to equilibrium. 

For the combined model, the virial radius (here defined as 
the radius where the mean enclosed density is 200 times the crit- 
ical density, assuming Hq = 71kms~ 1 Mpc~') is ^200 = 242kpc 

2 An animation and data comparisons for this ran 
can be found in the first author's presentation at 
www. ari .uni-heidelberg. de/meetings/milkyway2009/talks 
and at www. astro . umass . edu/~f ardal/movie_sphr_aug09 . gif . 
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Projected Radius (kpc) 

Figure 2. V-band surface brightness distribution of our model, com- 
pared to data. Solid magenta: total model surface brightess. Dot- 
dashed yellow: bulge. Triple-dot-dashed green: disk. Long-dashed pur- 
ple: halo. Dashed red: bulge and halo combined (or "spheroid"). Solid 
squares: Bro wn et al. (2003, 2007, 2 008) photometric fields. Open trian- 
gles: |prite^T&v^^nj3ergh ( 1994) photometric fields. Open diamonds: 
Gilbert et al. (in preparation) spectroscopic fields. All except the Gilbert 
points are taken from the minor axis. The inset expands the region studied 
in this paper; dashed vertical lines show the complete radial range of the 
target stars. 



and the corresponding virial mass is M 200 = 1.66 x 10 12 M Q . (For 
a virial threshold of 100, these quantities would be 3 19 kpc and 
1.91 x 10 12 Mq respectively.) The mass within 125 kpc, sometimes 
used as a benchmark of M31's halo mass, is 1.13 x 10 12 M Q . The 
mass of the stellar halo component within r2oo is 1.2 x 10 9 Mq, 
more than a factor of ten lower than the bulge mass. 

The surface brightness distribution produced by these com- 
ponents is shown in Fig. [2] Here we assume observed (not dust- 
corrected) M/Ly values of 5.6 for the bulge and 5. for the disk, 
based on estimates of their M/Lr and V — R color in lGeehan et al.l 
(2006). For the halo we assume a lower M/Ly = 2.5, which is rea- 
sonable for old stellar isochrones of the low metallicity ([Fe/H] < 
— 1.0) thought to be prevalent in the halo dKalirai et alJ 12006a: 
iKoch et ai] [2008). The smaller M/Ly of the halo compared to the 
bulge is justifiable based on the lower metallicity and extinction in 
the halo. Of course, these M/Ly values are fairly uncertain, and 
adopting a single M/Ly for each component is itself a gross over- 
simplification. 

There have been many observational studies of M31's sur- 
face brightness, most of which have concentrated on the SE mi- 
nor axis. We include three da tasets here: the photometric surv ey of 
Pritchet and van den Bergh jPritchet & van denB ergh 199^), the 
spectroscopic survey of Gilbert et al. (in preparation), and thre e 
points from deep HST surveys iBrown et al.l 00031 120071, 12008b . 
The Gilbert et al. values are calculated after removal of obvious 
substructure, and have been normalized to the surface brightness in 
the HST deep fields. The HST values in turn were calculated by 



summing the indi vidual stellar luminosities using data products in 
IBrown et ai] |2009l . This method should be more accurate because 
the HST fields probe nearly the entire luminosity function, because 
they require no uncertain sky subtraction, and because they directly 
overlap several of the spectroscopic fields. We calculate the normal- 
ization using three spectroscopic fields that overlap the HST fields, 
making a small correction for the slight differences in mean radii 
between corresponding HST and spectroscopic fields. We choose 
the normalization value using a minimum y} criterion, with results 
shown on the plot. It can be seen that the total surface brightness in 
our model traces observational data reasonably well, in particular 
reproducing the p r eviously known upwa rd kink around 20-30 kpc 
dlrwin et aljkoOSl ; ICourteau et ai]|201 ll) . Although some discrep- 
ancies may be present, these are no larger than the scatter in the 
observational points. Some of the observations are still likely con- 
taminated with substructure in the outer parts of the halo, so our 
model is chosen to follow the lower envelope of the data. 

Both bulge and halo components are represented here by 
spherical, isotropic populations. We regard the distinction between 
these hot bulge and halo components as somewhat artificial. The 
evidence for two components in the surface brightness data is 
merely the previously mentioned upward kink, and the separation 
between bulge and halo in our model is undoubtedly influenced by 
our choice of radial profiles. Thus we prefer to discuss them collec- 
tively as the "spheroid" component, which is shown by the dashed 
red line in Fig. [2] The transition between the surface mass density 
of these two components occurs at 25 kpc, and that in the V light at 
21 kpc. We note the decompositi on here is rather different than that 
for the /-band decomposition in ICourteau et al.l d201ll) . in which 
the bulge-halo transition occurs at a mere 3 kpc, due largely to dif- 
ferences in the assumed radial profiles. Courteau's preferred model 
uses a de Vaucouleurs bulge and a cored power-law halo, and these 
have very different dependences at small and large radii than our 
functional forms. The to tal spheroid compon ents in our model and 
in the preferred model of lCourteau et al.ld201 ll) are not so different, 
however. Even with a crude correction for V — I color, we find the 
two profiles differ usually by less than 0.3 magnitudes within the 
virial radius and never by as much as a magnitude. The two largest 
discrepancies come at the bulge-halo transition points in the two 
models, and lie either in the region where the disk dominates the 
spheroid, or in the halo region where the statistical and systematic 
uncertainties are high. 

The spheroid star velocity dispersion at 30 kpc in our model is 
l^knis" 1 and very close to Gaussian. From measurements else- 
where in M31's halo, we expect a nearly Gaussian distribution 
with a m ean of zero and dispersion of about 130k ms -1 in the hot 
spheroid dChapman et al. 2006; Gilber Tet al.l2007l) . suggesting the 
model spheroid component is reasonable. 

We use equal particle masses for all components of our simu- 
lation, so that the number density of simulation particles is propor- 
tional to the mass density within the model. In principle the ratio 
of mass density to the number of RGB stars that are selected in 
the observations is dependent on stellar population. However, the 
bulge and halo population at these radii have relatively old stellar 
popul ations lacking in stars with < 5 Gyr dBrown et al.l2003ll2006l 
2007), so their RGB density per unit mass should not be too dif- 
ferent from the stream material. As we will see later, the disk com- 
ponent plays at most a small role just in the innermost mask, so 
the possible large difference in its stellar population is not of great 
concern. Furthermore, accurate predictions of the kinematic com- 
ponent normalizations are never required by our analysis methods 
below. Rather, we allow freedom in fitting the normalizations of 
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Figure 3. Histogram of the stellar velocities in our sample. The solid his- 
togram shows only those stars classified as M31 giant stars, while the open 
histogram includes the contribution from foreground dwarfs. The dotted 
line shows the best Gaussian fit to the M31 stars assuming the expected 
spheroid velocity dispersion of 130 kms~', while the dashed line shows the 
best fit with velocity dispersion left free. Neither distribution adequately fits 
the data. 



the different components, and instead use their spatial and velocity 
distributions as a template to set the observed relative strength of 
the component contributions. 

When comparing to the observations, we enhance our parti- 
cle statistics by using a larger area than is covered by the observed 
masks. We choose particles within the region with major axis po- 
sition \X\ < 0.3°, slightly wider than any of our masks, as well as 
minor axis position 0.7° < Y < 2.2° . In this we rely on the fact that 
over small scales in the simulation the velocity distributions change 
slowly with position, and depend mainly on projected radius. 



3 ANALYSIS AND RESULTS 
3.1 Observed distributions 

In Fig.fJ] we show the overall radial velocity distribution of the ob- 
served stars. Throughout, we state velocity relative to M31's sys- 
temic velocity of — 300±4kms~ 1 dde V aucouleu rs~e~t al.ll99lh . We 
split the distribution into M31 and MW stars. As discussed in § 2, 
we expect a mean of zero and dispersion of ~ 130kms _1 if the 
distribution is dominated by the hot spheroid. Neither this distri- 
bution, nor the best Gaussian fit with dispersion left free, fits the 
distribution at all well, as can be seen in the figure and as veri- 
fied by % 2 tests. Most stars are contained within |v| < 50kms _1 , 
but a large fraction inhabit much broader tails of the distribution. 
In contrast, a fit with dispersion 130kms _1 to only the values with 
|v| > 100 kms -1 is statistically acceptable. The apparent peaks near 
— 175kms -1 and — 20kms~' hint at finer-grained structure. 

In Fig. [4] we show the distribution of our sample stars in 
metallicity. The stars are predominantly metal-rich, though a tail 
to low values is evident. In Fig. [5] we show the metallicity versus 
velocity of individual stars, color-coded by classification into M31 



Figure 4. Histogram of the stellar metallicities in our sample, with the solid 
histogram showing only the M31 giant stars. We use photometric metallic- 
ities throughout, since we believe these to be more reliable than spectro- 
scopic metallicities at the noise level of our spectra. 



vs MW stars. There is a small clump at — 175 km s which is the 
main contributor to the leftmost hump in Fig. [3] Note also that at 
velocities close to M31 systemic, there is a strong concentration 
of high-metallicity stars. This results in metallicity distributions of 
the stars with |v| < 75kms~ I and those with 75 < |v| < 150kms _1 
being inconsistent at the < 3% level, according to a Kolmogorov- 
Smirnov test. 

In Fig. [6] we show the metallicities of individual stars as a 
function of projected radius, with the median trend overlaid. The 
coverage gaps and variations in radial density are due to the posi- 
tioning and orientation of the masks: the two tangentially aligned 
masks produce narrow concentration of stars. The four radially 
aligned masks (with two overlapping) produce three broader but 
less dense zones of stars. The green line shows the giant-star metal- 
licities boxcar smoothed over their rank order in radius. This seems 
to indicate a rapid change in the metallicities around the edge of 
the shelf at about 1.83°, indicating a change in the mean popula- 
tion there. Comparing the metallicities of M3 1 stars in zones 4 and 
5, we find a 97% chance using bootstrap resampling that the mean 
is higher in zone 4, which supports the visual impression. 

In Fig. [7] we show the velocity versus radius of the sample 
stars. The uneven distribution in radius reflects the effect of the 
mask pattern. The velocity distribution of the M3 1 stars varies with 
radius. For example, a Kolmogorov-Smirnov test indicates that the 
distributions within radial zones 2 and 4 are inconsistent at the < 
10 -4 level. The velocity distribution is very far from Gaussian at 
most radii. Clumps of stars are obvious at about R = 1.33°, v = 
75kms~', andi? = 1. 77°, v = — 12kms _I . There also appears to be 
a wedge-shaped region within which many of the other stars reside, 
with these two clumps at two of its vertices. These velocity clumps 
are not especially concentrated in position within their masks, and 
are clearly not the result of an individual star cluster. The large 
number of M31 stars outside the wedge-shaped region also clearly 
shows a smoother hot spheroid is present. 
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Figure 5. Metallicity versus velocity for individual stars in our sample. Tri- 
angles denote M31 giant stars, and diamonds denote MW dwarf stars. (As- 
signed metallicities of the MW stars are not meaningful as explained in the 
text.) We have placed an ellipse around the small clump at v ~ — 175kms , 
[Fe/H] «s -0.75. 



Figure 6. Metallicity versus radius for individual stars in our sample, with 
symbols the same as Fig. \5\ The green line shows the metallicity boxcar 
smoothed over radial rank order. The numbers and separators near the top 
indicate the radial zones. 
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A small overdensity of 10 stars around R = 0.85°, v = 
89km s _1 is indicated by a green box. Observations bv lGeha et al.l 
fcOOfj) indicate NGC 205 has a tidal debris tail extending roughly 
in the direction of our innermost mask, and at their furthest de- 
tected extent of about 3.6 kpc from NGC 205's center the M31- 
centric velocity of this debris is about 80km s . The stars in the 
green box span a range of 3.2-3.8 kpc from NGC 205's center. The 
mean [Fe/H] of stars in the green box is —1.2. For comparison, the 
old stellar p opulation in NGC 205 has a mean of — 1 .06 ± 0.04 ac- 
cording to lButler & Martinez-Del gadol d2005l) . consistent with our 
stars within the uncertainties. Thus the kinematics and metallicity 
strongly suggest this overdensity is a continuation of a tidal tail 
from NGC 205. It remains to be seen whether the debris tail found 
here and in iGehaet all (2006) can be reconciled with the suggested 
kinematic detection of a t ail extending from NGC 205 to the NE 
dMcConnachie et alj|2004l) . 

This overdensity helps to explain the asymmetric velocity dis- 
tribution in our innermost mask, where the stars are mainly at pos- 
itive velocities. We note the innermost stars are at about 0.8° = 
48 kpc in the disk plane. It is possible that the disk contributes to 
some of the low-velocity stars in this mask, though a warp might 
be necessary to explain the offset from zero velocity. 

The strong wedge-shaped overdensity of stars is more signifi- 
cant and more challenging to explain in terms of known kinematic 
components of M31. Since the observations are taken along the 
minor axis, one would expect any contribution from M31's disk 
to have nearly zero veloc i ty. Ev en the extended cold disk compo- 
nent found bv llbataetal1d2005h . which has a substantial lag in ro- 
tation speed relative to the cold gas in the disk, cannot produce 
a 75kms~' line-of-sight velocity relative to M31, as exhibited at 
1.33°. The outer cold clump is at an appropriate velocity to be from 
the disk, but this location is at a projected radius of 24 kpc along 
the minor axis, or 107 kpc in M31's disk plane (for an inclination 
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Figure 7. Velocity versus radius for individual stars in our sample. Red tri- 
angles denote M31 giant stars, and cyan diamonds denote MW dwarf stars. 
The stars within the ellipse in Fig.|5]are circled. The green box indicates a 
clump of stars possibly associated with NGC 205. The blue polygon shows 
the test region used to calculate the velocity dispersion of the hot compo- 
nent. 
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of 77°). There is no sign of a disk contribution in the spectroscopic 
masks interior to the outer cold clump, or in the photometric survey 
maps (Fig. QJ, so we conclude the outer cold clump is not associ- 
ated with the disk either. 

As mentioned above, the disk contribution in the range 10- 
20 kpc on the SE minor axis is controversi al. Kinematic obser- 
vations suggest d ominance by the spheroid dKalirai et alj l2006a; 
Gilbert et al. 2007), while photometric observat ions have been in - 
terpreted as dominance by a disk component Jlbata et al1l2007l) . 
However, there seems to be agreement that the disk contribution 
is negligible beyond 20 kpc. Our kinematic results show this is also 
true from 18-24 kpc on the NW minor axis. 

To test the velocity dispersion of the hot component outside 
the triangular enhancement, we create a test region as shown in 
Fig. [7] We carry out a simple Bayesian analysis of the veloc- 
ity dispersion of the stars in this region, assuming they have a 
hot Gaussian velocity distribution centered at zero and with width 
<T V , and taking a flat prior for ln(o v ). This yields an estimate of 
o"V = 117± 15kms~', with the maximum likelihood value occur- 
ring at 113kms~'. Using a more conservative test region makes 
no significant change in these values. The hot component disper- 
sion found here is consistent with the dispersion a v ~ 130kms _1 
expected from previo us studies in other regions around M31 at 
slightly larger radii dChapman et al. 1 12006; Gilb ert et al.l 120071) . It 
is also in excellent agreement with the spheroidal component of 
M31 in our model, which has a dispersion ranging from 103 to 
118 km s _ 1 over the range 20-30 kpc . 

3.2 Comparison with model kinematics 

In Fig. [8] we again show the velocity versus radius of our stars, this 
time overlaid on the simulation points from the surrounding area. 
Purple points show the particles belonging to the GSS-related de- 
bris, and green the particles from the base M31 model. For slightly 
better agreement we have taken the liberty of reducing the radii of 
the GSS debris particles by 4%, an amount well within the model 
uncertainties (Fardal et al. in preparation). The agreement between 
the simulation prediction and observations is striking. The simu- 
lated GSS debris lies in a wedge that tapers to a tip at about 1 .9° or 
26 kpc, with a strong concentration at the upper edge. The wedge 
shape results from a caustic in projected ve locity that is a conse- 
quence of the kinematics of a radial shell (Merrifiel d & Kuiikenl 
1998). A large fraction of the observed stars are distributed in a 
similar pattern. 

In the rest of this subsection we will compare the distribution 
in the radius-velocity plane in more detail. Rather than just making 
a statistical comparison in arbitrary areas of this plane, we will try 
to interpret the observed results in a physical way. We first note that 
the slope of the observed wedge appears to agree with the simula- 
tion. S ince this slope is related to t he gravitational force at the shelf 
edge jMerrifield & Kuiikedfl998l) . this suggests our gravitational 
potential model near radius 20 kpc is reasonable. 

The wedge shape of the caustic derived in 
iMerrifield & Kuiiked dl998h is triangular, and symmetric around 
zero systemic velocity. However, this result assumes that the 
shell star orbits are exactly radial (no angular momentum) and 
monoenergetic (no velocity dispersion or energy gradient along 
the flow). In contrast, our simulated shelf material has both angular 
momentum and an energy gradient along the flow. In the Appendix 
we derive an analytic expression for the caustic envelope that 
incorporates both of these effects, as well as second-order terms 
in the potential and velocity projection. To use this expression we 
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Figure 8. Velocity versus radius as in Fig. UJ but this time including our 
simulation particles as small points. Particles are selected from a rectan- 
gle around the masks as discussed in § |2.2| Purple points show the shelf 
component formed by the GSS satellite debris, while green points show the 
component representing M3 1 itself, which is dominated by the hot smooth 
spheroid component. Red observed stars are to be compared to the combi- 
nation of green and purple model points. The orange lines show the analytic 
envelope discussed in the text. 

need parameters for the rotation and energy gradient which we 
now consider. 

The mean observed velocities are about — 12kms _1 in the 
large group of cold stars in zone 4 just inside the shelf tip, and 
— 18kms _I for the smaller group of six stars at slightly larger ra- 
dius and velocity near zero. Our simulation of the debris similarly 
assigns small negative velocities for these regions. We estimate the 
tip velocity in our simulation as about v t j p ~ — 20 km s~', by tak- 
ing a sample of particles at the shelf edge. In contrast, a shell model 
assuming purely radial motions for the stars would find a zero ve- 
locity relative to the host, at least to the measurement precision of 
the systemic velocity which here is roughly 4km s . The small 
negative velocity we find in the simulations has clear significance: 
it represents the effect of angular momentum within the shell parti- 
cles, which causes motion along the line of sight at the shell edge. 
We infer the same is true for the observed stars. The direction of the 
angular momentum results from the specific trajectory assigned to 
the GSS progenitor in the F07 model, in which stars after passing 
through the NE Shelf pass behind M31 into the W Shelf, curling 
toward us at the edge of the shelf. This trajectory is strongly con- 
strained by the need to place the stream and shelves correctly on 
the sky. 

In the simulation, the radial shell giving rise to the W Shelf 
can be seen to expand over time as stars of higher and higher orbital 
energies arrive there. This energy gradient is expected on quite gen- 
eral grounds, since stars with higher orbital energies take longer to 
complete their orbits. The energy gradient skews the wedge shape, 
increasing the apparent speed of outflowing stars (moving them fur- 
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Figure 9. Distribution of M31 stars in velocity within radial zones 2 (left), 3 (center), and 4 (right). Observations and model are kernel-smoothed with a 
Gaussian of a = 15kms~'. Black line: observed M31 stars. The shaded region indicates the lo~ error band. Magenta line: total distribution in model. Purple 
line: GSS-related model stars. Green line: M31 galaxy model stars. 



ther from M3 1 's velocity) and decreasing that of the inflowing stars. 
The simulated shelf stars mostly lie beyond the midplane, so that 
outflowing stars have positive velocities and inflowing ones have 
negative velocities; this results in an upward skew of the wedge. 

To combine the skewed velocity envelope given by equation 
IA4I with a global offset due to the angular momentum, we need 
four parameters. We estimate the shelf expansion speed v s /, = 
40km s _1 from the evolution of the shelf radius in the simula- 
tion. We use v rol = — 20kms~ 1 as the result of the line-of-sight 
projection of the angular momentum, as found above. We choose 
r s h — 1 .97° to match the observed three-dimensional shelf radius. 
Finally, from our potential model, we estimate the circular veloc- 
ity V c = [GM(r)/r] 1 / 2 to be 240 kms -1 at the shell radius. These 
choices yield the orange lines shown in Fig. [8] Clearly this enve- 
lope is a good match to the actual envelope of both the simulation 
and observed stars. In contrast, a symmetric linear wedge would be 
a poor fit to the simulated and observed caustic shapes. The upward 
skew of the wedge towards M31's center suggests that the simula- 
tion correctly places the bulk of the shelf beyond the midplane, so 
that the positive-velocity stars are the outflowing ones. 

In our simulation the M31 disk contributes a small amount 
to only the innermost mask, but in actuality it is not obvious that 
any stars with radii and velocities typical of the disk are present in 
our sample. Therefore our comparison to observations constrains 
only the shelf and spheroid components. We regard the particular 
bulge/halo decomposition of the spheroid in our model as some- 
what arbitrary, Formally, though, our M31 bulge component dom- 
inates the contribution to the surface mass density of the spheroid 
component out to 25kpc or 1.8°, in other words nearly the entire 
region probed by this survey, while the M3 1 halo component dom- 
inates outside that radius. 

In Fig. [9] and Fig. QJJ] we compare in greater detail the dis- 
tributions of simulated and observed stars in velocity space. We 
ignore zone 1 in all further comparisons with the simulation, due 
to the small number of detected stars and the apparent contamina- 
tion from NGC 205 debris and possibly M31's disk. (See Fig.[T]and 
Table[T]for the positions of the radial zones.) 

Fig.|9]shows smoothed-kernel representations of observed and 
simulated stars in the three central zones. In zone 2, there is clearly 
a steep-sided central component in both observations and simula- 
tions, weighted strongly to positive velocities. The simulation is 
more strongly peaked, and the smaller negative- velocity peak is 
not clearly distinct in the observations. This could be from an im- 
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Figure 10. Binned distribution of velocity difference Av from the caustic 
envelope shown in Fig. [8] as defined in the text. Points with error bars show 
the observed distribution. The magenta line shows the simulated distribu- 
tion. The significant jump at zero demonstrates the increase in density just 
inside the velocity caustic. 



balance between our spheroid and shelf components, or a slightly 
different spatial distribution of shelf stars in this range. In zone 3, 
the similarity between the simulation and the observations is re- 
markable. As we approach the edge of the shelf in zone 4, the cen- 
tral component shrinks to a cold peak. The stars in the observed 
and simulated cold components actually occupy a similar range, 
but a slightly different distribution of stars within this range pro- 
duces a small shift in the peak of the smoothed velocity distribu- 
tion. Clearly, the model provides an essentially correct description 
of the motions of stars within the shelf. 

The sloping caustics shown in Fig. [8] are smeared over veloc- 
ity in Fig.|9]by the finite range in radius within each zone. Fig. 1 101 
instead shows the velocity distribution in reference to this caus- 
tic velocity. We use a histogram rather than kernel smoothing to 
make a sharp division between stars inside and outside of the caus- 
tic. For each observed or simulated star i with 1.28° < Yj < 1.9° 
(i.e., from the inner edge of zone 2 to the shelf boundary), we com- 
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pute a quantity Av; = ±(v CTV ,i ~~ v i)> where v enV i ' s the nearer of the 
upper and lower caustics given by the orange lines in Fig. [8]at the 
star's radius. We use the plus (minus) sign for stars nearer the upper 
(lower) caustic. Thus Av,- is small and positive for stars just inside 
the caustic, and small and negative for stars just outside. We show 
the observed distribution as points with Poisson error bars, and also 
the result of the simulation as a solid histogram. It is clear that the 
density jumps by a large amount as one crosses the caustic bound- 
ary. The simulated distribution reproduces this distribution rather 
well, although the fit is poor from a statistical point of view. The 
same can be said of the Av distributions in individual fields which 
are not shown here. Not surprisingly, the shelf component is solely 
responsible for the jump at zero in the simulation, while the M31 
component has a broad, smooth distribution. 

Now that we have clearly established the presence of both 
hot smooth and shelf components from their velocity distributions, 
we use these distributions as templates to find the observed rela- 
tive strengths of these component. For each of our radial zones, 
we find the GSS and M31 velocity distributions at the velocities 
of the observed stars by kernel-smoothing the simulation veloci- 
ties. We combine these two distributions in proportion to their total 
masses within the zone, with the GSS component weighted by an 
additional factor g sat . We then compute the total likelihood lnL sat 
for the M31 stars being drawn from this distribution, by adding 
the logarithms of individual velocity distribution values. Repeating 
this process for a set of values for the weighting factor g sal , we 
find its maximum likelihood estimates and error estimates, and we 
can compare the likelihood lnL. sa , to that of the null model InLo, 
where we use only the M31 component and no GSS debris (equiv- 
alent to setting g sat = 0). This test offers very strong evidence for 
a satellite component in zones 2-4, whereas zones 1 and 5 are am- 
biguous. (This is not surprising since NGC 205 contaminates the 
innermost masks, while the outer mask should have very few shelf 
stars.) In our first set of tests we let g sat take on a different value 
for each zone. For zones 2, 3, and 4 respectively, we find optimum 
satellite weights lng sat of -0.8 ±0.5, -0.1 ±0.3, and 0.5 ±0.3. 
The log likelihoods relative to the null model are \n(L sar /Lo) of 8, 
29, and 31, indicating solid detections of the GSS component in 
each of these zones. If we assume a constant g sat and combine data 
from all fields, we find \ng sa i = 0.0 ±0.2, i.e. the optimum weight 
of the GSS component is g sa t — (1.0 ±0.2) the standard weight in 
our model. We note that ln(L Jflf /Lo) = 72, i.e. in a simple likeli- 
hood ratio interpretation the kinematic distribution is > 10 31 times 
more likely to be drawn from the model with a satellite component 
than that without. This overwhelming statistical requirement for a 
shelf kinematic component confirms the visual impressions from 
Figures |8land[T0l 

For future work on modeling the GSS, it will be useful to have 
an estimate of the contributio n of the shelf t o the M31 stars. The 
purely photometric sample of iTanaka et al.l d2010h shows a steep 
surface brightness drop at the W Shelf edge, which we can esti- 
mate to be roughly 3 magnitudes in their most metal-rich channel 
by extrapolating their outer surface brightness profile inwards. This 
would suggest the spheroid constitutes approximately 0.06 of the 
total metal-rich stars and the shelf component itself about 0.94, al- 
beit with large uncertainties. 

3 One should not read too much into the fact g sa , is near 1 here, as we 
used the considerable observational freedom in the spheroid modeling of 
§ 2 to choose a model where this was assured. If we normalize the GSS 
component by this weighting factor, however, all of our other results are 
virtually independent of which spheroid model we use. 



The ratio of shelf to spheroid components varies with position, 
and the average ratio obtained with our somewhat irregular spatial 
coverage is of no particular interest. Therefore we use the simu- 
lations to estimate this ratio averaged over an area defined by the 
radial range 1.35-1.95° and the azumuthal range 0-30° to the SE 
of the NW minor axis. We compute the satellite-contributed frac- 
tion f sat from the numbers of simulation particles in this region in 
the GSS and M31 components, weighting the GSS count by the 
satellite weight g sal . Using the likelihood as a function of g sat just 
described, we find the GSS stars make up f sat = 0.59 ± 0.05 of 
the total stars in this region. Better definition of the GSS-related 
structures can be obtained by focusing on the metal-rich popula- 
tion. If we repeat the likelihood calculation using only metal-rich 
stars with [Fe/H] > —0.75 to constrain g sat , we require a higher 
weight of g S ai = 1.34 ±0.37. This yields a higher satellite fraction 
of fsat = 0.84 ± 0.05 for this subset of stars. (This suggests that 
the shelf is more metal-rich than the spheroid, which we discuss in 
more detail later.) These estimates are the most precise determina- 
tions yet of the ratio between the shelf and spheroid, and should be 
very useful in constraining future models of the GSS. 

The observed velocities appear somewhat asymmetric about 
the origin; most noticeably, the strong enhancement at R = 1.33°, 
v = 75kms~' is not matched by an equally strong clump at 
— 75kms~ 1 . The model shelf component shows an asymmetry in 
the same direction (compare the strengths of the upper and lower 
caustics in Fi gure s [8l and 1 1 0) . The simulated shelf stars lie beyond 
the midplane, and we have argued above the observed stars must as 
well due to the skew of the wedge in r-v space. Therefore the pos- 
itive velocity stars are outgoing and the negative velocity stars are 
infalling. The ratio of these two components tests the density gra- 
dient along the stream, since the infalling stars have larger orbital 
phase. 

To test the model versus the data in a more precise way we 
restrict our attention to zones 2 and 3, which have a large fraction 
of shelf stars but are far enough away from the edge that the precise 
dividing line between outgoing and infalling stars is not important. 
In the simulation, the outgoing stars (classed as such by their actual 
3-d radial velocities) make up (67 ± 1)% of the shelf stars within 
this region. The dividing line between outgoing and infalling stars 
is reached here at — 8kms _1 . We use the GSS and M31 velocity 
distributions, along with the maximum likelihood satellite debris 
weight factor of f sw = 0.7 1 , to assign a GSS mixing fraction to each 
observed star. We can then estimate the number of outgoing and 
infalling GSS stars by summing up these mixing fractions above 
and below the dividing line of — 8kms _1 . The observed outgoing 
fraction is 70 ± 3%, in excellent agreement with the simulation. The 
density therefore falls off steeply with increasing orbital phase. If 
we accept that the stream density should rise when approaching the 
phase of the progenitor, then this demonstrates the shelf lies at a 
larger orbital phase than the main body of the progenitor, as in our 
model. 



3.3 Metallicity of shelf and spheroid components 

Our maximum-likelihood weighting of GSS and M31 components 
over all zones produces for each observed star a mixing fraction, 
that is, a probability that the star is from one component or the 
other. We now use these mixing fractions to estimate the metallicity 
distributions of the individual components. 

In Fig. [TT] we show the metallicity distribution for the entire 
sample of M31 stars, kernel-smoothed with a Gaussian of disper- 
sion 0.1 dex. Dashed colored lines show the contribution of the 
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Figure 11. Kernel-smoothed distribution of our W Shelf sample stars in 
metallicity (solid black curve). The dashed curves indicate the contribution 
of the shelf (purple) and M3 1 proper (green) components, using the initial 
decomposition based on kinematics. The solid curves with the same colors 
indicate the contributions using the mixing fractions based on both kine- 
matics and metallicity. The observed metallicity distribution of stars within 
the GSS itself, abo ut 100 kpc to the so uth, is shown by the orange kernel- 
smoothed curve forlGilbertetal.l f2009). and by the squares with error bars 
for llbataetai]j2007ll 



individual components, where each star contributes an amount pro- 
portional to its mixing fraction for a given component and is kernel- 
smoothed as before. It is clear that the GSS debris has higher metal- 
licity on average than the smooth spheroid. 

Our kinematic separation of the stellar components is far from 
complete — there are many stars with significant mixing fractions 
in both components, meaning the metallicity distributions of the 
components are still to some degree mixed together. This mixing 
can be reversed by applying the metallicity distribution as an addi- 
tional criterion in computing the mixing fractions. Here we assume 
that metallicity is independent of position and velocity throughout 
our observed regionQ Starting with the original mixing fractions 
from kinematics alone, we compute the metallicity distributions of 
each component and use it as an additional factor in the likelihood 
of a star being from that component. We then recompute kernel- 
smoothed metallicity distributions from the new mixing fractions, 
and iterate this process until the mixing fractions converge. In 
technical jargon, this procedure involves a non-parametric mix- 
ture model with prior fuzzy classifications on all the sample points, 
which we solve for the mixture components using an expectation- 
maximization algorithm. The result is shown as the solid purple 
and green curves in Fig. [TJJ This algorithm always amplifies the 
differences between the original component distributions. It natu- 
rally amplifies noise as well, so we constructed mock samples from 
our simulation of similar size to the observational sample and tested 

4 A constant metallicity distribution in the observed phase space is prob- 
ably a good assumption for the W Shelf material itself, since the stars are 
moving rapidly through the radial shell. For the M31 model, the range of 
radii we sample is less than a factor of 2 and the spheroid dominates over 
the entire range, so the distribution in this component is probably reason- 
ably constant as well. 



the separation using the remaining simulation particles to define the 
components. We found that at our sample size, the algorithm does 
a good job of recovering the original component distributions with- 
out introducing much visible noise. 

Our final metallicity distributions in Fig. [TT] show an even 
clearer separation of the satellite debris and smooth spheroid com- 
ponents. The mean metallicities are —0.7 for the satellite debris 
and —1.2 for the spheroid component, for a typical difference of 
0.5 dex in [Fe/H]. The spheroid metallicity has a similar mean and 
dispersion to that found in the smooth, dynamically hot spheroid 
population at slightly larger rad ii (30-33 kpc) on the opposite side 
of M31 bv lKaliriri et alH2006j) . 

The figure also illustrates the distribution o f metallicity 
within the GSS itself, from the kinematic sample of Gilbert et al.l 
J2007I) and the gro und-based photometry of llbata et alj d2007h . The 
Gilbe rtet al.l d2007h sample has advantages of kinematic selection 
and identical reduction procedures to those employed in our W 
Shelf sample, but in truth the two GSS distributions are quite sim- 
ilar. The W Shelf debris component shows remarkable agreement 
with these GSS samples, obtained in a region ~ 100 kpc away. 



4 DISCUSSION 

The observations presented here reveal two main components be- 
sides the foreground contaminants from the MW: a hot smooth 
spheroid with a roughly Gaussian velocity profile, and a shelf com- 
ponent with kinematics matching that of a radial shell. The kine- 
matic properties of the shelf component match very well to a sim- 
ulation of the GSS's formation. The metallicity distributions of the 
two components are distinct, with the shelf having higher mean 
metallicity and matching very well to observations of the GSS. 
This agrees with the earlier conclusion of lTanaka~e t al. (2010j), who 
compared purely photometric samples in the W Shelf and GSS re- 
gions and also found them to be nearly identical. 

Some caution is warranted when drawing connections be- 
tween tidal structures based on their stellar populations. Different 
tidal progenitors can have similar metallicity distributions, and a 
single tidal structure can have gradients along it, as seems to be the 
case for the Sagittarius stream dChou et al .1120071) . Indeed, the GSS 
itself seems to have a metallicit y gradient from it s core to its less 
dense "cocoon" on the SE side (I bata et al.ll2007|) . The models of 
iFardal et alj d2008l) reproduce such a gradient by assuming a gra- 
dient in the progenitor. Importantly, however, neither observations 
nor models suggest a strong gradient along the stream, especially 
for regions separated by a comparable amount from the progeni- 
tor's core. In the F07 model, the W Shelf leads and the GSS trails 
the current core location by similar amounts. Thus it is reasonable 
to connect the two structures by their stellar populations. 

The precision of resolved-star spectroscopy allows a remark- 
able range of conclusions to be drawn from the velocity structure. 
We showed the observations agree not just with the overall wedge 
pattern of the shelf component, but also agree with more subtle as- 
pects. The ratio of upper and lower caustic populations, the velocity 
shift at the tip of the wedge from M31's systemic velocity, and the 
slightly skewed shape of the wedge all have physical meaning in the 
simulations, and the data are in good agreement with the simulation 
trends for all these aspects. The high precision of the kinematic data 
encourages use in constraining our satellite model further. For ex- 
ample, we already had to adjust the radii in the model very slightly, 
an adjustment well within the constraints posed by other observa- 
tions of the stream debris. 
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Together, the metallicity and kinematics make a strong case 
that the W Shelf indeed originates from the same disruption event 
that produced the GSS. The structure around M31 is still incom- 
pletely understood. In particular, part of the GSS region shows a 
second stream offset by + 100km s from the GSS itself, which 
is currently unexplained by our models. Is it possible that the W 
Shelf itself results from a similar, but independent accretion event 
to that which formed the GSS? The rich variety of observational 
features matched by our simulations suggests that this is unlikely. 
However, further confirmation of the model is worthwhile. The next 
test for the model is likely to come in the NE Shelf region, where 
the model places the bulk of the GSS-related material. In this re- 
gion our model predicts a similar wedge kinematic pattern to that 
in the W Shelf (see F07). We have argued that negative-velocity 
(outflowing) stars have been detected already in this region, but the 
detection of a positive-velocity (returning) component will provide 
a clear connection between all parts of the GSS model. The progen- 
itor's total stellar mass in our model is roughly comparable to that 
of the LMC, and adding up the stellar mass of all the GSS-related 
features will provide another test of the model. 

The smooth spheroid component is consistent with that found 
elsewhere in M31. The velocity distribution agrees in its mean and 
dispersion with that found previously on the opposite side of M31's 
disk, as does the metallicity distribution. This provides strong ev- 
idence for a symmetric, mostly virialized hot component at inter- 
mediate radii. It may be a bit simplistic to call this inner spheroid 
"smooth", as evidenced by the — 175kms~' clump of stars with a 
small range of metallicities we found in §[3] This is not unexpected 
in theories of galaxy halo formation, as fine-grained structure can 
linger in phase space and stel lar population space far longer than is 
appar ent in morphology (e.g. lHelmi & Whitelll999t I Johnston et al.l 
2008). However, at a minimum it is made up of a very large num- 
ber of components below the threshold of individual detectability 
spread over a large velocity range, as opposed to the cold mono- 
lithic kinematic pattern of the shelf component. 

It is clear that M31's inner halo remains a confusing region, 
and much work remains to be done to decompose it into kinemat- 
ically distinct components. The extent of the W Shelf in particular 
remains to be mapped, as in star-count maps it becomes confused 
with M31's outer disk on both the NW and SE sides. Extrapolating 
from progress in the last decade, we expect rapid growth in our un- 
derstanding of M31's halo from the combination of large imaging 
surveys with spectroscopic surveys such as we have presented here. 

The wedge pattern detected in Fig. [8] is to our know ledge 
the fir st time a radial shell kinematic pattern dMerrifield & Kuiikenl 
has been clearly detected in any galaxy, although features 
of the pattern were suggested earlier in F07. It is also the first 
spiral galaxy found to exhibit a system of multiple shell struc- 
tures, though recently NGC 4651 has been found to exhibit sim- 
ilar structures (Martmez-Delga do et al . 2010). Shell systems are of 
course common in elliptical galaxies (e.g. | ]Vlalin & Carted [l983l ; 
iMerrifield & Kuiikerj|l998l : ICanalizo et al.il2007n . Our results offer 
inspiration for the future detection of shell kinematics in more dis- 
tant galaxies and its use in measuring their gravitational potential, 
using either bright tracers such as planetary nebulae or globulars or 
else RGB stars detected with a future generation of telescopes. 
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APPENDIX A: PHYSICS OF LESS-IDEALIZED SHELLS 
— EXPANSION, ROTATION, AND VELOCITY 
DISPERSION 

Here we derive some results on the kinematics of non-ideal shells, 
for use both in the preceding analysis of the W Shelf and in under- 
stan ding future observations of ot her stellar shells. 

iMerrifield & Kui jken ( 1998) derived the observable kinemat- 
ics of a stellar shell in the idealized case of a monoenergetic, per- 
fectly radial flow of stars in a gravitational field of uniform strength. 
In the same approximation, F07 derived the shell surface density in 
terms of the gravitational potential and flow rate of stars through 
the shell. Both papers discussed several physical and observational 
aspects of stellar shells besides those mentioned in this paper. How- 
ever, here we need several results not mentioned in either paper. 

To establish n otatio n, we repeat the argument of 
IMerrifield & Kuiiked H998t) here. Denote the three-dimensional 
radius as r and the projected radius as R. Consider a shell of 
fixed radius r s / t , with stars of constant orbital energy flowing 
out to r s f, and back, and for convenience let us consider only 
the stars on the far side of the host. Write b = r/r s i„ /3 = R/r s / r 
Let the gravitational field have a constant strength g; then at 
the shell edge the circular velocity is V c = (gr s h) 1 ' 2 - The radial 
velocity of the stars is given by v 2 = 2g(r s f, — r) = 2V 2 {\ — b). 
The line-of-sight velocity is instead given by v/ os = (z/r)v r , 
where z is the distance relative to the host distance. (Here we 
neglect perspective effects.) Near the shell edge we can to 
first order take z/r = [(r 2 - R 2 )/r 2 ] l l 2 w [2(& — /3)] 1/2 . Then 
v /o! = ±2V c (b - /3) 1/2 (1 - b) x l 2 . The magnitude of the velocity 
is maximized at b = (1 + /3)/2, implying the envelope of v; OJ is 
given by v env = ±V C (1 — fi) = ±V c (r s ft — R)/r s h. Here the top and 
bottom signs apply to the outgoing and incoming stars respectively. 
Thus projecting the velocities has squashed the sideways parabola 
of v r (Y) into a wedge. The gradual turnover in v; 0J around its 
maximum gives the distribution of stars in velocity space a caustic 
or hom shape. A similar argument applies to stars on the near side 
of the host, with the roles of the positive and negative caustics 
reversed. 

While undoubtedly useful as a first approximation, this con- 
tains one fundamental problem: the radial shell exists only because 
stars have different orbital energies and hence arrive at the shell 
edge at different times. The gradient is invariably in the sense that 
later-arriving stars have higher energies, because larger orbital en- 
ergy increases the orbital period. We can handle this problem sim- 
ply by incrementing the outward velocity of all stars in the shell by 
the same amount. 

(vr-vsh? =2g(r sh -r) = 2V 2 {l-b) (Al) 

If we label each star by the time t p it arrives at the shell boundary, 
its specific energy becomes E = v 2 h /2 + + v^gtp. The energy 
gradient can then be expressed as 

dE/d{t p ) = v sh g = v sh V 2 / r sh . (A2) 

(As argued above, dE/d(t p ) and v s / t are both in general positive.) 
This shows the expansion speed of the shell is proportional to the 
gradient in orbital energy along the stream. 

We can again solve for the maximum line-of-sight velocity 
at a given projected radius. Defining e = v s i 1 /{gr s f l ) x l 2 = v^/Vc, 
S = 1 -/3, and v = (e/8)([165 + e 2 ] l l 2 - e), and differentiating to 
find the extremum, we find an envelope velocity of 

v W = V c (5±v) 1 / 2 [e±(5 T v) 1 / 2 ]. (A3) 

This results in an envelope in r-v/ OJ space that is shifted upwards 
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relative to the monoenergetic case at most radii, with an increas- 
ingly skewed shape approaching the shelf radius so that the tip 
remains at zero velocity. This is a natural result of the outgo- 
ing stars (upper branch of the envelope) having higher energies 
and thus greater velocities. For lines of sight outside the radius 
where the "returning" particles have zero radial velocity, given by 
(25) I//2 = e, there is no returning caustic and the minimum enve- 
lope velocity of zero is reached at the midplane. The tip remains 
fixed at zero velocity because the radial velocity at the shelf edge 
lies perpendicular to the line of sight. If the shell material instead 
lies on the near side of the host galaxy, then the energy gradient 
instead skews the envelope in the negative-velocity direction, with 
the tip again fixed at zero. If the shell material fills both sides of the 
galaxy, the degeneracy of the near and far side envelopes is broken 
and both envelopes will be present, offset from each other at all 
r<r sh . 

Since the expansion velocity is typically small compared to 
V c , it is more useful to take a perturbative approach, which allows 
us to incorporate second-order effects in the gravitational potential 
and velocity projection. To this order we can write zjr ~ V2(S — 
£>) 1//2 (1 + \D- ^5). with D = 1 -b. We assume a flat rotation 
curve, in which case <t>(r) - <P{r sh ) w - W}D{ 1 + \D) . To this order 
we can assume the extremum is reached at D = S /2, the lowest- 
order solution. Combining the various factors, we reach our final 
approximation 

(A4) 

The first factor of (1 + h8) is contributed by the projection factor, 
and the second (in the case of a flat rotation curve only) by the 
potential. 

In the case of the W Shelf, as will be true in most cases, the 
shell material is not purely radial but instead has significant angu- 
lar momentum. This affects the observed velocities by shifting the 
envelope up or down. We apply this as a global shift to the entire 
envelope in Fig. [8] A higher-order analysis would take into account 
the radial dependence of tangential velocity and derive the enve- 
lope shift by rederiving the caustic analysis above, but the simple 
shift used here appears accurate enough for our purposes. 

As a final exercise, we estimate the smearing in the stellar ve- 
locities near the caustic edge induced by a local radial velocity dis- 
persion a r of the stream particles. This dispersion might be left over 
from initial conditions in a young stream like that discus sed here, 
or els e result from interaction with LCDM substructure dCarlberd 
2009) in an older one. To first order we can estimate the observed 
dispersion from the monoenergetic picture, where the projection 
factor at the velocity caustic is z/r = (1 — j3) 1//2 . Thus we expect 
the envelope edge to be blurred by an amount 

Olo S = {{rsh-R)lr sh )\ llZ <yr- (A5) 



V c 8 l l 2 {\ + -8) e±8 x l 2 {\ + -8) 



© 2012 RAS, MNRAS 000,|JJ{T4| 



